Generalized coherent state representation of Bose-Einstein condensates 
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We show that the quantum many-body state of Bose-Einstein condensates(BEC) consistent 
with the time-dependent Hartree-Fock-Bogoliubov (TDHFB) equations is a generalized coherent 
state (GCS). At zero temperature, the non-condensate density and the anomalous non-condensate 
correlation are not independent, allowing us to eliminate one of the three variables in the TDHFB. 



I. INTRODUCTION 



The recent experimental realization of Bose-Einstein condensates (BEC) in supercooled trapped atoms has stimu- 
lated great interest in the theoretical description of the quantum state of BEC [1-6]. The exact quantum state and 
, many body wavefunction of this system is not known, and simple approximations such as the Hartree approximation 
are commonly employed. Knowing the quantum state is crucial for describing the BEC dynamics; in particular, the 
I ' many-body hierarchy which leads to an infinite sequence of progressively higher order equations may be truncated 
*^ consistently by making an assumption on the quantum state of the system [7] . 

Despite the lack of an exact representation of the quantum state, it has been possible to make significant progress 
by using physically sound assumptions about how certain operator products should be factorized in order to truncate 
the many-body hierarchy. The Gross-Pitaevskii Equation (GPE) has been used for describing the dynamics of zero 
temperature trapped, atomic BEC with great success [8-11]. The time-dependent Hartree-Fock-Bogoliubov (TDHFB) 
equations [12-15] are coupled nonlinear equations connecting the dynamics of the condensate and non-condensate 
atoms required for the description of finite temperature BEC in the collisionless regime. Other finite temperature 
• theories include the time-dependent Bogoliubov-de Gennes equations [16], the Quantum Kinetic Theory [17-22], and 
' Stochastic methods [23-27]. However none of these treatments directly addresses the precise quantum state of BEC 
f— i . that consists of the condensate as well as the non-condensate atoms. We note that the condensate atoms are described 
at the mean field level in various theories. However, the GPE is derived under the assumption that all correlations of 
annihilation and creation (c|) operators for the non-condensate atoms in some basis state i vanish e.g.: 



{c\c)c k c m ) = (cjcfc) = (c]c m ) = (1) 



o 

S while the HFB equations are derived using a different ansatz: 

' {c\c]c k c m ) = {c\c k ) (c]c m ) + (c^c k )(c\c m ) + (cjcj) (c k c m ) . (2) 

o 

It is straightforward to show that the GPE may be derived by assuming that the quantum state of BEC's at zero 
temperature is a coherent state [2]. This is closely related to the description of the laser by a coherent state in 
Quantum Optics [28,29]. Indeed one of the earliest stated goals in BEC research has been the development of an 
"atom laser," the matter- wave equivalent of laser [30-32]. It has been argued that, owing to the presence of the 
intrinsic interatomic collisions, the zero temperature BEC is more accurately represented by a squeezed state rather 
than a coherent state [5,6]. 

Eq. (2) used for deriving the HFB equations is reminiscent of Wick's theorem for a system in thermal equilibrium 
[7,33-37]. This implies that the thermal equilibrium state described by a statistical density matrix is clearly a possible 
candidate for the quantum state of BEC. However, this choice does not provide a satisfactory physical picture for 
the TDHFB equations that describe dynamical condensates away from equilibrium. In addition, collision-induced 
squeezing [5] has not been included in the current description of finite temperature BEC. The identification of a 
quantum state that can describe the dynamics of finite temperature BEC, consistent with the TDHFB equations is 
thus an open issue. 

In this paper, we propose a generalized coherent state (GCS) ansatz for the many body density matrix describing 
the dynamical quantum state of BEC. Such states were originally used to describe anharmonic dynamical systems 
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such as many body interacting fcrmions/bosons [38] while preserving some of the useful properties of the original 
Glauber's coherent states for the harmonic oscillator [28,29]. They encompass the Glauber coherent state as well as 
the squeezed state as special cases. The GCS are particularly convenient for formulating variational dynamics because 
of certain algebraic structures originating from the underlying Lie group algebra [38,39]. Using this ansatz, we derive 
the TDHFB equations [13] via the time-dependent variational principle. This principle which allows the description 
of the many-body system in terms of a small number of parameters is intimately related to classical Hamiltonian 
Poisson bracket mechanics that describe classical dynamics from the minimum action principle. 

The paper is organized as follows: In Section II, we review the key properties of GCS relevant to variational 
dynamics at zero and finite temperature. In Section III we derive finite-temperature variational equations of motion 
in both the real space and the trap basis. Conclusions and discussion are given in Section IV. 



II. GENERALIZED COHERENT STATE VARIATIONAL DYNAMICS 



Mathematically, a set of GCS is determined by a Lie group G, its irreducible unitary vector representation T with 
the space V and a reference state |0) € V. The GCS are defined as states that have a form T(g)\Q) with jeG. 
More specifically, for a generic quadratic Hamiltonian in some operators Tj 

H = ^c i f i + Y^c ij f i f j , (3) 

i i,j 

the Lie group G is characterized by the commutation relations amongst the complete set of operators Tf. 

[f i ,f j } = Y / C* j f k , (4) 
k 

where are known as the structure constants of the set {Ti}. 

For a harmonic oscillator, {Ti} — {di,d],i} where hi, a] are the boson annihilation and creation operators and / 

is the identity operator, and C*- = giving the ordinary Heisenberg-Weyl group. On the other hand, the operator 

set {Ti} = {a\aj,a\a}j,aia,j,i} may be used to construct the following Hamiltonian that describes the system of 
many-body interacting bosons: 

H= ^2 Hijal&j + ^2 Vijkia\ & j & kai- (5) 

ij ijkl 

An extended Heisenberg-Weyl algebra may be obtained by a repeated application of the standard boson commutators 
[Sj,a]] = Sij. Writing T-^ = hihj, = djd], = a\a,j + the non-vanishing commutation relations that 

define the extended Heisenberg-Weyl algebra are [39] : 

\f{-) f (+)i - x fW i a f(z) + x tW+a f( z ) (61 

\f (z) f , (-)l — _A f (-) _ A T(-) (7) 
Kmn) J rs J — U mr- L ns U ms- L nr \ < J 

Pmni ^rs~^] = S nr T^J + S ns T^J , (8) 

while the mixed commutators between the linear and bilinear operators have the form: 

Pmn > Qj] = — finj&tn (9) 

[f£l& j ] = -6 mj a n (10) 

[^rnn i — ^nj^m H~ ^rnj^n (H) 

[f^,a]} = 6 nj al. (12) 

At zero temperature, the unnormalized generalized coherent states that belong to such extended Heisenberg-Weyl 
algebra for interacting bosons have the form [38]: 

|V(t)> =exp I ^>;(r)a] +E^( T K ta ] ) l fi o>- (13) 
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The reference state \£lo), which can be normalized to unity (f2o|^o) = 1, may be chosen arbitrarily. However, 
construction of a useful set of coherent states for a given dynamic system depends crucially on the choice of |17o) 
which also determines the structure of the phase space of the dynamical system [38]. We shall take Eq. (13) to be our 
generalized coherent state ansatz for BEC where the reference state |fio) is the particle vacuum state of the ordinary 
Heisenberg-Weyl group i.e. di\fl ) = 0. The goal of the variational approach is to determine the time-dependent 
parameters aj(r) and Pij(r) that represent the evolution of the state |^>(t)) . Eq. (13) indicates that the GCS has 
a form which combines a coherent state and a squeezed state. A coherent state is given by applying a displacement 
operator D(a) on the vacuum \a) = £>(a)|0) where D(a) = exp(^ i a*ai — cticn). Since exp (— a*a,i) |0) = |0), an 

unnormalizcd coherent state may be written as |a) = exp |0)- O n the other hand, a unitary squeezing 

operator is given by — exp ^ CijO-i^j ~ £ij a ! a ]) > where ^ = exp (z%) is an arbitrary complex number. 

GCS zero temperature variational dynamics is obtained by implementing the dynamical variational principle as- 
suming that the space of trial wavefunctions M is represented by a set of GCS. One possible way of formulating the 
variational dynamics in Hilbcrt space is based on projecting the vector Hx for any x G M (H being the Hamiltonian 
operator) into the tangent subspace to M at x. This leads to a vector field in M that determines the variational 
dynamics. 

The variational equations at zero temperature are derived as follows: Given a Hamiltonian H, and time-dependent 
wave functions |^(r)), we minimize the action: 



S[Q(t)] = J dr\i (n(T)\dn(T)/dT) - (Q(t)\H\Q(t))^ . (14) 



By choosing a GCS form for |fi(r)), the resulting variational equations can be written in the Hamiltonian form for 
any set % of coordinates which parametrize 

^ = (15) 
where {• • •} denote Poisson brackets and H is the classical Hamiltonian defined by: 

H(Q) = (Sl\H\Q). (16) 

The use of Poisson brackets clearly establishes the link between the variational equations and the classical dynamics. 
When the classical Hamiltonian is given by 

« = EE M: ) .. in <4)---mj, a?) 

n— 1 i\---i n 

the Poisson bracket assumes a very simple form provided the wave functions \Q) are parametrized by the expectation 
values (f2|Tj|f2) of the operators Tj rather than by the parameters fij. These expectation values then constitute a 
full set of parameters that uniquely specify the quantum state In particular, if the operators Tj form a closed 
algebra Eq. (4) , the Poisson brackets for Tj is given by: 

{T m ,f n }=i^2 C ^nTk, (18) 

k 

and the variational equations of motion for T m take the closed form 

^ = EE E cu^AU^-'-ifa d9) 

n—l j — 1 io---i n 

For the Hamiltonian Eq. (5), it therefore suffices to calculate the equations of motion for the expectation values (a\a,j) 
and (aiCLj) to uniquely specify the dynamics. In the derivation of the equations of motion we use the differential 
property of the Poisson brackets: 

{/, 9h} - ~{gh, /} = {/, g}h + g{f, h}. (20) 

When the expectation values are for generators of the set of GCS of some Lie group G, their Poisson brackets are 
given by the commutators of the underlying generators of the group. This simplifies the calculation greatly, as it 
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gives direct correspondence between the ordinary quantum mechanical commutators and the Poisson brackets. The 
variational procedure is then formally equivalent to the Heisenberg equations of motion. It should be noted that the 
transformation between the expectation values and the parameters may be tedious. However the transformation is 
never used explicitly; suffice it to know that such transformation exists and we can then proceed to derive closed 
equations of motion for the expectation values. 

The formulation of variational dynamics at finite temperatures constitutes a more complicated task for the following 
three reasons: (i) At finite temperatures the system evolves in Liouville space and a set of trial density matrices Ml 
rather than wavefunctions needs to be identified, (ii) an attempt to project Lp with p e Ml and L is the Liouville 
operator into the tangent space to our ansatz faces a difficulty since the Liouville space does not have a natural 
scalar product that can be used for this projection (the scalar product inherited from Hilbert space i.e. the overlap 
of two density matrices does not have a direct physical significance), (iii) The compatibility of the equilibrium and 
the dynamical approaches that is straightforward at zero temperature (i.e. the state in Ml with the lowest energy 
must be the stationary point of the dynamical equation) is not so obvious in Liouville space. In Appendix A we show 
that all of these issues can be adequately resolved for GCS. In particular we show that, the trial density matrices 
may be assumed to have the form of finite temperature equilibrium density mar ices, and that the equilibrium and 
the dynamical approaches are compatible so that the trial density matrix which minimizes the Helmholtz free energy 
is, indeed, a stationary point of the dynamical equations. These results will be used in the coming section. 



III. VARIATIONAL EQUATIONS FOR INTERACTING BOSONS 



A. GCS in real space 



In this section we apply our formalism to a system of interacting bosons described by the following Hamiltonian 
with a fixed chemical potential p: 



H = Ho+H f (t) 



(21) 



with 



/ 



Ho = / dr^(r) 



1 

2^ 



A + V t rap(r) - p 



Mr) 



+ drdr'ft(r)ft(r')V(r-r')i>(r')i>(r) 
H f (t) - / dv^{v)V s (v)^{v). 



(22) 
(23) 

where Vt rap (r, t) is the magnetic potential that confines the atoms and Vf (r, t) denotes a general time- and position- 
dependent external driving potential. An infinite-dimensional extended Heisenberg- Weyl algebra is generated by the 
operators V'(r), ip\ Y(r,r') = -0(r)^(r'), yt(r,r'), N(r,r') = ^ (r)tp(r'), and /. The space of the representation 
in which the Hamiltonian [Eq. (22)] is defined can be described as the space of wave-functionals ^[a;(r)], where the 
one-particle boson operators are given by: 



5x(r) 



x(r) 



Sx(r) 



x(r) 



(24) 



x(r) is a harmonic oscillator coordinate associated with position r and Eq. (24) can be used to represent of the 
operators Y, Y\ and N. 

At zero temperature the set of coherent states represented by Gaussian wavefunctions 



V[x(r)] = Aexpj-i J drdr'a(r,r')[x(r) - x (r)][x(r') - x (r')] 



(25) 



are parametrized by complex-valued functions xo(r), and cr(r,r'). According to the formalism developed in Section 
II and Appendix A the coherent finite temperature density matrices p[xj(r)] with j = L,R [Left (ket), Right (bra)] 
are represented by Gaussian wavepackets: 



p[xj(r)]=Z ^xpj-i^ J drdr'a kj (r,r')[x k (r) 



xf(r)][x 3 (r')-xf\r')} 



(26) 
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These trial density matrices are parametrized by the vector and matrix functions x^\r) and (Tkj(r, r'). We note 
that the Gaussian wavepacket [Eq. (26)] constitutes a coordinate representation for a density matrix of the form 
p = Z^ 1 cxp(-K), where K is given by a combination of linear and bilinear terms in single-particle operators. 

The variational parameters of Eqs. (25) and (26) which denote the displacement and the width of the Gaussian 
wave packet in phase space are related to the average number of particles and the quantum mechanical squeezing of 
the number-phase conjugate variables in BEC. These parameters may be related to physical quantities such as the 
condensate fraction and the excitation energy, by transforming Eq. (25) or Eq. (26) to the quasiparticle basis; the 
resulting relationship between the variational parameters in different bases is not simple. However this transformation 
is never used explicitly since the GCS ansatz allows us to derive equations of motion directly for the parameters of 
interest; the expectation values of the relevant operators. 

The most convenient parametrization for trial wavefunctions [Eq. (25)] or density matrices [Eq. (26)] is given by 
the expectation values of linear and bilinear combinations of boson single particle operators: 

z(r) = Mr)), «(r, r') = (Y(r, r')) - z(r)z(r'), p(r, r') = (N(r, r')> - z*(r)z(r') (27) 

where the expectation value is taken with respect to the wavefunctions given by Eq. (25) or density matrices of Eq. 
(26). Wick's theorem [7,33-37] allows us to express the expectation value of any operator in terms of the parameters 
given in Eq. (27) both at zero and finite temperature. The dynamical equations for the system of interacting bosons 
may therefore be derived in the same way for both zero and non-zero temperatures by starting with the Heisenberg 
equations of motions for linear and bilinear combinations of the single-particle operators and then evaluating the right 
hand sides using the Wick's theorem. This results in closed equations for the parameters z(r), n(r, r'), and p(r, r'), 
which will be derived next. 

B. Variational equations of motion in real space 

The Heisenberg equation of motion for ip(r) reads: 



dt 

where 



H sp $(r) + J dr'^(r')V(r,r')^(r')tP(r) (28) 



if sp (r) = -^A + V t ra P (r), (29) 

nr,r')E^[V(r-r') + y(r'-r)], (30) 
and we have used the commutation relations for the boson field operators 

$(r),^(r / )]=J(r-r / ), [^(r)^(r')] = [^( r ), fitf)] = 0. (31) 

Taking the expectation values of Eq. (28) and noting the definition of /c(r, r') and p(r,r') [Eq. (27)], we obtain the 
equation of motion for the mean field: 

l h^->-=H sp {r)z{r) + J dr'V{r,r') {\z{r')\ 2 z(r) + z*{r') K (r,r') 

+ z(v')p(r,r') + z(v)p(v',v')} + V f (r,t)z(r). (32) 

The equations of motion for n(r, r') and p(r,r') can be derived similarly by computing the time derivatives using 
Eq. (27) and Eq. (28) in the product rule: 

lh ^dP~ = HSP ^ r ^ + J d *"V(r'y) {i(v",v')p(v,v")+i(v",v")p(v,v') 

+ C>',r>*(r",r)} - H sp {v')p{v, r') - J dr"V(r,r") {f(r, r")p(r", r') 
+ |(r", r")p(r, r') + C*(r, r>(r", r')} + V f (r, t)p(r, r') 

-V f (r',t)p(r,r') (33) 
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ck(r,r') 

in ; 

dt 



= /T*(r)«(r, r') + J dr"V(r', r") {|(r", r>(r, r") + f(r", r>(r, r') 
+ C(r',r")[ i o*(r,r") + <5(r-r")]}+i/ sp (r') K (r,r')+ ^ dr"^(r,r") 
x {|(r",r) K (r',r") + |(r",r")«(r,r') + C(r,r")p(r",r')} 



+ V f (r,t)K(r,r') + V f (v',t) K (v,v'), 

where we have introduced the auxiliary functions 

£(r,r')=z*(r)z(r') + p(r,r') : 
C(r,r')=z(r)z(r') + K (r,r'). 



(34) 

(35) 
(36) 



For the commonly used special case of the contact interatomic interaction for V(r, r'), Eqs. (32-34) are simplified 
greatly; these are given in Appendix B. 

An important consequence of the GCS ansatz is that at zero temperature the functions p(r, r') and k(t, r') arc, in 
fact, not independent [40]. By deriving an explicit relationship between them, it is possible to eliminate the p(r, r') 
variables. This relation is derived in the trap basis in Appendix C, and then converted into the real space basis: 



P(r,r') 



f>(r,r')+ / K*(r, y")k(t" , r')dr" — —S(r, r') 



(37) 



It can be verified by direct substitution that once Eq. (37) holds initially, it remains true throughout the dynamical 
evolution. The reduced set of equations are then the coupled equations Eqs. (32) and (34) with p(r, r') replaced by 
the expression Eq. (37). The two independent variables z(r) and k(t, r') constitute a very convenient parametrization 
of squeezed states. z(r) represents the average position whereas n(r,r') are responsible for squeezing. This can be 
easily understood from the fact that the coherent state is an eigenstatc of the annihilation operator ip(r) while a 
squeezed state is generated using a squeezing operator which is a function of the quadratic operator Y(r, r') in the 
extended Heisenberg-Weyl algebra. 



C. Variational equations of motion in the trap basis 



For completeness, we outline below the derivation of the same variational equations in the trap basis. In this basis, 
the Hamiltonian is written as: 



H = Hijal&j + Y v Hkia\a]a k ai + ^ E H a \°-j- 



ij ijkl ij 

The matrix elements of the single particle Hamiltonian Hij are given by 



Ha = |rf 3 r0*(r) 



-<£A + V trap( r) 



</>j( r ), 



(38) 



(39) 



where the basis state (j>i(r) is arbitrary; a convenient basis for trapped BEC is the eigenstates of the trap since is 
then diagonal. The indices may also be viewed as the mode indices in a multimode quantum state. The symmetrized 
two particle interaction matrix elements are 



Viju = 7: 



(ij\V\kl) + (ji\V\kl) 



where 



(ij\V\kl) = yd 3 rdV0*(r)0*(r')nr-r')^(r')0/(r), 
with V(r — r') being a general interatomic potential. Also, 

En= J rfr^(r)V>(r,t)^(r), 



(40) 



(41) 



(42) 
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where V/(r, t) denotes a general time- and position-dependent external driving potential as denned previously. First, 
we proceed by establishing that our generalized coherent state Eq. (13) at zero temperature is a Gaussian in coordinate 
space. With the choice of the particle vacuum state of ordinary Heisenberg-Weyl algebra as our reference state \Qo), 
Eq. (13) the action of the operators <Xj and a\ on the wave function Cl(qi, . . .qjy) in the coordinate representation, 
where Af is the total number of bosons is: 



«i = 



a = 



~V2 V% 



(43) 



and the conditions aifl(qi, . . . q^) = for i = 1 . . . Af imply that the reference state £lo(qi, ■ ■ ■ qu) is a Gaussian in 
coordinate space: 



^o(gi, ■ •• qx) = 



exp 



air) 



(44) 



It should be noted that, consistent with our choice of basis in Eq. (38), the index i of the coordinate variable qi in 
Eqs. (43) and (44) refers to the trap basis "mode" i. In addition, although we are using the trap basis, the finite 
total number of particles Af implies that the vector space used is effectively a finite dimensional space spanned by a 
truncated set of trap basis states. Acting on this wave function with the generalized displacement operator[Eq. (13)] 
preserves its Gaussian form since the action by the operators hi and a\ simply shifts the origin while the operators 
ciitij and a\a\ change the variance. The resulting state is thus a Gaussian of the form: 



&{q\ , ■ ■ ■ qu) = A exp 



(45) 



where r/j , i = 1 . . . Af are the complex numbers which determine the average position while oij is an Af x Af symmetric 
matrix that determines the covariances or the amount of squeezing. For an ordinary coherent state, ov,- 
Gaussian with unit covariance. 

Similarly, at finite temperatures, the trial density matrix takes the form (k,j = L,R): 



Sij i.e. a 



I kj Im 



sf (0) ][<4-d 0) i 



(46) 



The GCS ansatz may therefore be considered to be a ground state of some effective quadratic Hamiltonian in 
coordinate and momentum operators. For such a Hamiltonian, any correlation function can be represented in a path- 
integral form where the action only has linear and bilinear terms [41]. The resulting Wick's theorem is then identical 
to that for a thermal state. 

The expectation values to be used in the parametrization our state are the condensate mean field z i; the non- 
condensate density Pij and the non-condensate correlations : 



Zl = (Si) p i:j = (ataj) - {a]){aj) = {a^aj) - {ai){aj). 



'At 



(47) 



Since these variables are the expectation values of the generators of the set of generalized coherent states of the 
extended Heisenberg-Weyl algebra, their Poisson brackets are given by the commutators of the underlying generators. 
This results in TDHFB equations of motion for Zi , and Kij . 

At zero temperature, the relationship between p^ and [Eq. (37)] is: 



'1 



1 



(48) 



This enables us to reduce the number of equations. More details of this relation are provided in Appendix C while 
the TDHFB equations in the trap basis including the simplified zero temperature form are given in Appendix D. 
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IV. DISCUSSION 



Using the GCS ansatz, we have derived variationally the TDHFB equations equations of motion for BEC, which 
are known to be valid in the collisionlcss regime. This implies that the GCS ansatz should be applicable in the lower 
temperature, collisionlcss regime. It should be noted that the HFB theory has several inconsistencies such as the 
violation of the Hugenholtz-Pines theorem [42] which states that the excitation spectrum should be gapless in the 
homogeneous limit [12]. This issue has been addressed by various authors; for instance, the Popov approximation, in 
which the anomalous correlation is neglected, was shown to give a gapless spectrum [12]. Recently, it has been shown 
that by replacing the contact interaction potential with a more sophisticated pseudopotcntial, many of the inconsis- 
tency problems of the HFB equations including the violation of the Hugenholtz-Pines theorem, inconsistencies with 
the many body T-matrix calculations, and the ultraviolet ! divergences can be overcome [43]. In this paper we have 
presented our HFB equations with the general interaction in both the real space and the trap basis; pseudopotentials 
such as those discussed in Rcf. [43] can thus be accommodated. 

Since the GCS is a squeezed state, the present work may be considered an extension of a previous result that 
demonstrated stationary BEC to be squeezed [5] and a more recent result that has shown that dynamically evolving 
BEC under the time-dependent GPE described using the Hartree approximation (i.e. pure condensate, no non- 
condensate atoms) is squeezed [6]. 

The representation of the dynamical quantum state of BEC as a GCS provides physical insight about the total 
system of condensates plus non-condensates in terms of particle annihilation and creation operators, and how it 
evolves as a whole in the Schrodinger picture. The interdependence of p and k at zero temperature enables us to 
eliminate the p variables from these equations, reducing the size of the problem and simplifying the numerical solution. 

The multimode squeezing, and hence the entangled state nature of BEC is clear from the form of the quantum 
state, Eq. (13). The study of quantum entanglement is currently gaining great interest owing to its importance 
in quantum information theory [44,45] as well as in the understanding of the foundations of quantum mechanics 
[46,47]. Experimental sources of generalized coherent states in matter- waves already exist in the form of atomic 
BEC's. However, ways to access and manipulate this type of matter-wave entanglement remains an open challenge. 
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APPENDIX A: PROPERTIES OF FINITE TEMPERATURE GCS 

As noted in the main text, the formulation of variational dynamics at finite temperatures requires us to address the 
following issues: (i) At finite temperatures the system evolves in Liouville space and a set of trial density matrices Ml 
need to be identified, (ii) an attempt to project Lp with p € Ml and L is the Liouville operator into the tangent space 
to our ansatz faces a difficulty since the Liouville space does not have a natural scalar product that can be used for 
this projection, (iii) The compatibility of the equilibrium and the dynamical approaches is not so obvious in Liouville 
space. In this Appendix, we show how all these issues may be adequately addressed. We start by introducing the 
trial density matrices. 

Let A be the real Lie algebra of the real Lie group G involved in the definition of a set of GCS. For our case, A 
is the basic (real) algebra generated by the generators cij, a], and T mn ; the elements a belonging to the algebra A, 

a e A, are then linear combinations of these generators. In addition, let A^ and be the complexification of A 
and G i.e. G^ is the complex Lie group that corresponds to the complex Lie algebra A^ c \ Complexification of an 
algebra (group) gives an algebra (group) generated by the original generators for which the coefficients are allowed to 
be complex, rather than real numbers. More specifically, it means constructing a complex analytical algebra (group), 
for which the "real" version is the original one. For example, given a real structure which is a map p : G — > G the 
real part of the group consists of points g so that p(g) — g. If one defines p(g) — g^ one has G rea i — SU(2), while if 
one defines p(g) = g* , one has G rea i = SL(2, R). SL(2, C) is then a 3-dimensional complex analytical group (3 is its 
complex dimension) which serves as complexification for both SU(2) and SL(2,R). 

The representation T of a group G associates with any g £ G a linear operator T(g) (acting in some complex vector 
space referred to as the space of the representation) so that T(g 2 gi) = T(g 2 )T(gi). T can be naturally extended to a 
representation of G^ in the same vector space V. A representation of a group has the corresponding representation 
of an algebra and vice versa, and in the corresponding algebra representation we have T{a) being operators in the 
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same space for a e A with T(ai +a 2 ) = T(ai) + T(a 2 ) and T([ai, a 2 ]) = T(ai)T(a2) — T(a 2 )T(ai). The representation 
of ^4 can also be easily extended to representation of A^ . 

We define the manifold Ml of normalized trial density matrices represented by p{g) = Z^ 1 {g)T{g) for all elements 
g belonging to the complexification of G, g € G^ c \ so that T(g) is hermitian where Z(g) = TrT(g). The normalization 
condition Trp(g) = 1 is obviously satisfied. In the GCS case the projection that closes the dynamical equation can 
be formulated as follows: We define a tangent vector v(p) that satisfies the following property: 

Tr{T(a)v(p)} = Tr{T(a)Lp} (Al) 

for all a e A^ and p e Ml- For most practical applications there is one and only one tangent vector v for any p that 
satisfies Eq. (Al). In this case, a well-defined vector field v(p) describes the variational dynamics in the manifold Ml, 
of trial density matrices. The physical meaning of Eq. (Al) is clear: If we refer to the operators T(a) with a e A^ 
as the fundamental operators, the variational dynamics is obtained by the requirement that the dynamical equations 
hold for the expectation values of the fundamental operators. This implies that, similar to the zero temperature case, 
the variational equation of motion at finite temperature is derived using the Heisenberg equations of motion. 

We assume that the trial density matrices are represented by finite temperature equilibrium density matrices with 
Hamiltonians given by the fundamental operators: 

p = Z- 1 exp[-(3T(a)] = Z~ l T{g) (A2) 

for some a G A^ c \ This implies g = exp(/3a) i.e. elements g € G of the Lie group G can be represented as the 
exponentials of the corresponding complex Lie algebra elements a and it follows then that g G G^ c \ because if A is 
the Lie algebra of the group G, then A^ is the Lie algebra that corresponds to G^ c \ 

We conclude this section by demonstrating that this way of closing the dynamical equation [Eq. (Al)] guarantees 
the compatibility of the equilibrium and dynamical variational approaches. Using the variational approach, the 
equilibrium density matrix can be obtained by finding the minimum of the Hclmholtz free energy 

F(p) = Tr(Hp) - p-'Trip log(p)) (A3) 

among the normalized trial density matrices p e Ml- The requirement SF = yields: 

SF(p) = Tr(HSp) - P^TriSpp) (A4) 

for any tangent 8 p. It follows from Eq. (Al) that Sp = [a, p] is tangent for any a e A^ c \ This yields: 

Tr{T(a)L(p )} - Tr{T(a)[H, p Q }} - Tr{H[T(a), p }} 

= /3- 1 Tr{logpo[T(a),p ]} = /3- 1 Tr{T(a)[logp , Po}} = (A5) 

which implies that v(po) = 0. Stated differently, the trial density matrix which minimizes the free energy is a 
stationary point of the dynamical equations. 



APPENDIX B: VARIATIONAL EQUATIONS IN REAL SPACE FOR THE CONTACT POTENTIAL 



For contact interatomic interaction, V(r, r') = UaS(r — r'), Uo — 47r ^ a , where a is the s-wave scattering length and 
m is the mass of a single atom. This implies that the integrations in Eqs. (32-34) are removed: 



ih^- = H s P(r)z(r) + U {\z(r)\ 2 z(v) + 2p(v, r)z(r) + k(t, r)z*(r)} + V f (r, t)z(r) 



ih 



dp(r,r') 
dt 



dt 



H s P(r) +2tVf(r',r')J p(r, r') + U ((r' , t')k*(t, r') 

H s P(r>) + 2U i(v, r)] p(v, r') - U (*(r, r)«(r, r') 
-V f (r,t)p(ry)-V f (r',t)p(ry) 

H s P(r) + 2U i(v', r')] k(t, r') + Uofc , v')p*(v, r') + U C(r', r') 

+ [H s P(r') + 2U £*(t, r)] «(r, r') + U ((r, r)p(r, r') 
+ V f (r,t) K (r,r') + V f (r',t)K(r,r'), 



(Bl) 



(B2) 



(B3) 
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where £(r, r) and C(r,r) are as given in Eqs. (35-36). The equations of motion for p(r, r',t) and k(r,r',t) may be 
written in the compact form: 

ih^- = Eg - QYt (B4) 

where we have defined 2x2 matrices 

T,( T1 >\-( ^ r ' r ') A(r',r') \ C (rr')-f /?(r ' r ' ) K(r ' r,) ^ fB^ 

E(r ' r) -i-A*(r,r) -fc*(r,r')J ^ ' } ~ I «*(r, r') p*(r,r') + lj' (Bd) 

and 

h(r,r') = H s *(r) + V f (r,t) + 2U £(r' , r'), (B6) 

A(r,r) = (7oC(r,r). (B7) 

Eqs. (Bl) and (B4) constitute the TDHFB equations for the contact interatomic potential approximation, in real 
space [12]. 

APPENDIX C: ZERO TEMPERATURE RELATIONSHIP BETWEEN p AND k 

In order to derive the relation between pij and Ky at zero temperature using the GCS ansatz, we find that it suffices 
to consider GCS state such that (0|aj|0) = i.e. Gaussian wave functions centered at q = 0. These states form 
an orbit M of the group G which corresponds to the algebra generated by Tmn and T$nn ■ We shall introduce a set of 
functions Smn, Smn on M 



SW(fi) = (fi|fW|fi), fiW(fi) = (fi|fi±)|fi) (CI) 
and define two sets of auxiliary functions 

-Pmn(^) = ^ [^ma ^c«r^ — ^mL^an ) (^2) 

a 

G mn (fi) = 2 [si^W - SWSW] . (C3) 

a 

Our aim is to show that F mn (Q) is a constant i.e. its derivatives are zero. In particular, showing that F mn (Q) — 5 mn 

and identifying the expectation values Smn and Smn, in terms of p^ and Kij for our example completes the required 
proof. It is found that the derivatives of F mn (Q.) are linear combinations of G mn (f2), and therefore it suffices to prove 
that the auxiliary function G mn (f2) is zero for all m and n. 

We note that since [fyKf^] = 0, the operators f^ +) which are considered as vector fields on M determine a 
complex structure on M. A function / is said to be holomorphic if it satisfies the condition Tmn f = 0. Operator 
Tmn then represent derivatives in the antiholomorphic direction. 

In particular, the functions Smn constitute a set of holomorphic coordinates in the vicinity of £1 where Qo represent 
states with Gaussian wave functions. It can be shown that Smn(Qo) — while Smn(^o) — S mn so that 

^mn(^o) = &mn- (C4) 

A direct calculation yields T^G mn — which implies that G m n is holomorphic and can therefore be written as a 

series in S^ in the vicinity of the point S^ = (i.e. £1 ). Since Sml(Q) = S ma , it follows from Eq. (C3) that the 
expansion of G mn starts with the second-order terms: 

oc 

Gmn=J2 G rnl ( C5 ) 

We next define the degree of a function / by Df — degf = \ . T 1 ^ f . It is clear that degS^ = ±1, degSff = 0, 
and deg(fg) == degf + degg. It follows from Eq. (C3) that degG TOn = 1. On the other hand Eq. (C5) implies that 
dcgGmn = j and therefore contains the degrees of 2 and higher. This implies that G mn = for all m and n. 



10 



It can be verified that T^m F mn is a linear combination of G a b and hence TmnF mn = 0. Similarly, by conjugating 
the relation Tmn F mn = 0, Tmn F mn = 0, implying that F mn (Q) is a constant. This, together with Eqs. (C2) and 
(C4) imply 



^mn\, il -) — / J \ °ma°an a ma°a 



(C6) 



Since S ( m 2 = K mn and Sml = p mn + \5 mn Eq. (C6) gi 



ivcs 



+ 2^') 2 e K i P K Pj — 



(C7) 



Solving for p^ finally yields: 



Pij 



rSij + 



(C8) 



or in matrix form, 



where J is the unit operator. 

Eq. (C9) may also be expanded in a Taylor series as: 



p = k^k — (k^k) 2 + 4 (k^k) 3 — 



(C9) 



(CIO) 



It is possible to gain further insight into the nature of GCS from the fact that Eq. (C9) holds if the quantum state 
of the system is a quasiparticle vacuum state \0) qp such that 0i(t)\O) qp = where $i(t) is the quasiparticle annihilation 

operator of the Bogoliubov transformation, cn(t) — J2j^o Ujij3j{t) + Vj*/3j(i). For a quasiparticle vacuum state, 
and Kij may be written in terms of matrices U and V as follows: 



Pij — 2-*l "pi "PI 

p/0 



V r *Vni and 



p#0 



which implies 



p 2 +p = n^n and pn — np*, 



(Cll) 



(C12) 



using the orthogonality and symmetry conditions between the matrices U and V, UU^ — VV^ = 1 and UV T — VU T = 0. 
These relations can be shown to be identical to Eq. (C9) by solving the quadratic equation in p. The quasiparticle 
vacuum state may therefore be considered a squeezed state of condensate and non-condensate atoms. 



APPENDIX D: VARIATIONAL EQUATIONS IN THE TRAP BASIS 

1. TDHFB Equations 



The TDHFB equations in trap basis is given as follows: 

ih—TT = E Hij Z J + E ^i kl [ z j z k z l + ZpjkZl] + ^2 VijklKklZ* + ^2 EijZj 



dt 
dt 



jki 



kl 



E 

r 

+E 



H ir + 2 ^2 Viklr {z%zi + pik) + E Vl 



kl 



Pro - E 



H r j + 2 ^2 Vrjki {zlzi + pik) + E r 



kl 



r L kl 



^ Virkl {zkZl + Kkl) 



K rj E 



^2v rjk i {zlzf + K* kl ) 



kl 



(Dl) 



(D2) 
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dt £-~t 



H ir + 2 Viklr {Zk Z l + Plk) + Ei- 



kl 



K rj + 



H r j + 2 ^2 Vrjkl {zkz'i + p* lk ) + E r j 



kl 



E 



^2V irk i (z k Zi + K k l) 



kl 



Prj + Yl 



^ Vrjkl (z k Zi + K k l) 



Pir 

+ ^2V ijk l (z k Zt+K kl ). 



(D3) 



kl 



The TDHFB equations in real space derived in Section III, Eqs. (32-34), may be transformed to the corresponding 
trap basis, Eqs. (D1-D3) in a straightforward manner, using the following relations between the real space basis and 
the trap basis variables: 



z(r) =^2zi<j>i(r), p{r,r') = ^ p ij( j>* (r)<^ (r'), k{y,y') = ^ «y ;4>i{*)<t>j (V), 



(D4) 



along with the definition of the tetradic matrix V^ki 

1 



Vijki = 



(ij\V\M) + (ji\V\M) 



where 



(ij\V\M) = Jd 3 rd 3 r' ^(r)0*(r')^(r-r')0 fe (r')0 ; (i 



(D5) 



(D6) 



with V(r — r') being a general interatomic potential. Under the contact interaction approximation, Vijki takes a 
simpler form: 



Vij k i — 



m 



dr#(r)#(r)^(r)Mr). 



(D7) 



2. TDHFB at zero temperature 



The TDHFB equations, Eqs. (Dl) - (D3) hold for all temperatures. However, we have noted that the variables p 
and k are not independent variables for the generalized coherent state ansatz at T = 0. pij can therefore be eliminated 
using the following relation 



-I + kU- -I 

4 2 



+ o~P, 



where the function Sp = for T = 0, the TDHFB equations take the form: 
dz 



where 



ih-rr = H z z + H„z* +Ez 
at 



1 , 1 T . 

-1 + K'K 1 + Op 

4 2 



ih^- = [h, + k^k - ^7 + Sp] - (kA* - An 



A + A 

ih- 



l r ( 1 - 

—I + K'K 1 + dp 

4 2 



+ A, 



1 



dn 



-K + K 1 — 

dt 



[Hzjtj = H i3 + ^2 Vkij 
ki 

[Hz*]ij = VjkijKki 



z t z i + 2, jSih + Y K tm K ™k - 5i k + 2Spi k 



ki 



hij = Hij + 2 J2 V 



kij 



kl 



z k z l + \ J S lk + ^2 K lm K ™k - ^Sl k +Spik 



Ei. 



Aij = ^2 Vijki [zkZi + Kki] ■ 



ki 



(D8) 

(D9) 
(D10) 

(Dll) 

(D12) 
(D13) 

(D14) 

(D15) 
(D16) 
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At T = 0, Sp = and p = jl + k^k — |7. By direct differentiation of p 



1 



/'//^ = '7' r _ - , 

dt VT+Ik^k \ dt 



dn 



K + K dt 



while Eq. (D2) implies 



// t ^ = [/ M -/ --/]-(/vA -A/v ). 



(D17) 



(D18) 



Since Eqs.(D17) and (D18) are equivalent, both the left and the right hand side of Eq. (Dll) are zero at T = 0, 
5/7 = i.e. the only independent equations to be solved are Eqs. (D9) and (D10). 

For comparison, we note that the GPE, which is a zero temperature theory for a coherent state ansatz, is simply 
obtained from Eq. (D9) by setting pij — Hij = 0: 



dt 4^ 



kl 



(D19) 
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